;;########################################################################
;; unimob2.lsp
;; contains code to implement univariate analysis model object
;; file contains interpration, report, visualize and analysis methods
;; Copyright (c) 1995-8 by Forrest W. Young
;;########################################################################

;3 ViSta system methods (continued)


(defmeth univar-model-object-proto :interpret-model ()
  (let* ((name (send (send self :data-object) :name))
         (df    (send self :t-df))
         (t-val (send self :t-value))
         (p-lev (send self :t-p-level))
         (p-lev-string)
         (n-ways (send self :n-ways))
         (eta-sq (/ (^ t-val 2) (+ (^ t-val 2) df)))
         (real-resp-var "the sample mean")
         (null-phrase "the sample mean")
         (meandif (mean (send self :var)))
         (real-group-var)
         (direction (send self :direction))
         (ci-level  (send self :ci-level))
         (ci-limits (send self :t-ci-limits))
         (mu (send self :mu))
         (anal-type "single sample")
         (mean-str "mean")
         (compare "one") (nulsign "equals") (sign "does not equal") 
         (zt (if (send self :sigma) "Z" "T"))
         )
    (cond 
      ((send self :paired) 
       (setf mean-str "mean difference")
       (setf meandif  (mean (send self :var1)))
       (setf anal-type "paired (dependent) samples")
       (setf real-resp-var "the differences between the paired observations")
       (setf null-phrase "the difference in sample means")
       )
      ((> n-ways 0) 
       (setf mean-str "mean difference")
       (setf meandif (- (mean (send self :var1)) (mean (send self :var2))))
       (setf anal-type "independent samples")
       (setf real-resp-var (first (send self :real-varnames)))
       (setf null-phrase "the difference in sample means")
       (setf real-group-var (second (send self :real-varnames)))
       ))
    (cond
      ((= direction 0)
       (setf compare "two"))
      ((> direction 0)
       (setf nulsign "is less than or equal to")
       (setf sign "is greater than"))
      (t 
       (setf nulsign "is greater than or equal to")
       (setf sign "is less than")))
    (if (< p-lev .0001) 
        (setf p-lev-string "p < .0001") 
        (setf p-lev-string (format nil "p = ~6,4f" (fuzz p-lev 4))))
    (setf str (format nil "Interpreting ~a univariate tests~2%" anal-type))

    (unless (send self :sigma)        
            (setf str (strcat str (format nil "CONFIDENCE INTERVAL:~%For the ~a data we are ~2d% confident that the interval from ~4,2f to ~4,2f includes the true population value. For other samples (all the same size as for these data) we would expect the confidence interval to randomly move back and forth such that 95% of all possible samples would have a 95% confidence interval that would include the true population value.~2%" name (truncate (* 100 ci-level)) (first ci-limits) (second ci-limits) ))))
    (setf str (strcat str (format nil "~a-TEST:~%" zt)))
    (setf str (strcat str (format nil "For the ~a data, the ~a-tailed, ~a ~a-Test with ~d degrees of freedom, yields ~a=~3,2f with ~6,4f. This means that when the null hypothesis that ~a ~a ~3,2f is true, the probability is ~6,4f of obtaining a ~a as extreme or more extreme than the ~a of ~6,2f observed here.~2%It also means that over many samples (all the same size as for these data), the probability is ~6,4f that we would be wrong when we reject the null hypothesis that ~a ~a ~3,2f in favor of the alternative hypothesis that ~a ~a ~3,2f.~2%" name
                                  compare anal-type  zt df zt 
                                  (fuzz t-val 2) p-lev-string
                                  null-phrase nulsign mu  
                                  p-lev-string 
                                  mean-str mean-str meandif
                                  p-lev-string  null-phrase nulsign mu
                                  null-phrase sign mu)))
    (unless (or (equal "mean" mean-str) (send self :sigma))
            (when (or (send self :paired) (not (> n-ways 1)))
                  (setf str (strcat str (format nil "STRENGTH OF RELATIONSHIP: ~%For the ~a data we estimate that ~2d% of the variability in ~a is related to sample membership. As with the other statistics interpreted above, this estimate will vary from sample to sample." 
               name (truncate (* 100 (fuzz eta-sq 2))) real-resp-var)))))
    (help-window str :title "Interpreting Univariate Tests")
    ))

(defmeth univar-model-object-proto :report (&key (dialog nil))
  (let* ((var (send self :var))
         (var1 (send self :var1))
         (var2 (send self :var2))
         (var1-label (send self :var1-label))
         (var2-label (send self :var2-label))
         (real-varnames (send self :real-varnames))
         (real-varnames-string)
         (real-respvarname)
         (n-ways (send self :n-ways))
         (mu  (send self :mu))
         (direction (send self :direction))
         (ci-level  (send self :ci-level))
         (ci-limits (send self :t-ci-limits))
         (df    (send self :t-df))
         (t-val (send self :t-value))
         (p-lev (send self :t-p-level))
         (p-lev-string)
         (n1 (length var1))
         (mean1  (mean var1))
         (stdev1 (standard-deviation var1))
         (n2 nil)
         (mean2 nil)
         (stdev2 nil)
         (meandif nil)
         (stdevdif nil)
         (sterrdif nil)
         (pooled-variance nil)
         (dflist nil)
         (compare "one")
         (sign "NE  ")
         (nulsign "EQ  ")
         (paired-st "A SINGLE SAMPLE")
         (effect-string nil)
         (zt "T")
         (sigma (send self :sigma))
         (w nil)
         (mus "Mu")
         (nonpar-st "WILCOXON'S SIGNED RANK TEST FOR")
         (nonpar-au "Wilcoxon's   T")
         (F-test nil) 
         (F-df nil) 
         (F-p nil)
         )
    
    (when var2
          (setf n2 (length var2))
          (setf mean2  (mean var2))
          (setf stdev2 (standard-deviation var2)))
    (cond
      ((= direction 0)
       (setf compare "two"))
      ((> direction 0)
       (setf nulsign "LTEQ")
       (setf sign "GT  "))
      (t 
       (setf nulsign "GTEQ")
       (setf sign "LT  ")))
    (when (send self :sigma) (setf zt "Z"))
    (when (send self :paired) 
          (setf meandif  (mean var))
          (setf stdevdif (standard-deviation var))
          (setf sterrdif (/ stdevdif (sqrt (1+ df))))
          (setf paired-st "TWO PAIRED (DEPENDENT) SAMPLES"))
    (when (> n-ways 0) 
          (setf paired-st "TWO INDEPENDENT SAMPLES")
          (setf nonpar-st "MANN-WHITNEY TEST FOR")
          (setf nonpar-au "Mann-Whitney U")
          (setf mus "Mu1 - Mu2")
          (setf sterrdif (send self :std-err-diff))
          (setf pooled-variance (send self :pooled-variance))
          )
    (setf w (report-header (send self :title)))
    (display-string (format nil "UNIVARIATE ANALYSIS of the ~a Data."
                            (send (send self :data-object) :name)) w)
    (setf real-varnames-string
          (if (= (length real-varnames) 2) 
              (strcat "s " (first real-varnames) " and " (second real-varnames))
              (strcat " "  (first real-varnames))))
    (display-string (format nil "~%Analysis based on variable~a"
                           real-varnames-string) w)
                               
    (display-string (format nil "~2%~a-TEST FOR ~a:" zt paired-st) w)
    (display-string 
     (format nil "~%Null hypothesis:        ~a ~a ~7,3f" mus nulsign mu) w)
    (display-string
     (format nil "~%Alternative hypothesis: ~a ~a ~7,3f (~a tailed)" 
             mus sign mu compare) w)
    (when sigma (display-string 
       (format nil "~%Known Population Standard Deviation:   ~7,3f" sigma) w))
    (display-string (format nil "~2%SAMPLE STATISTICS:") w)
    (when (> (length var1-label) 12) 
          (setf var1-label (subseq var1-label 0 12)))
    (when (> (length var2-label) 12) 
          (setf var2-label (subseq var2-label 0 12)))
    (setf var1-label (subseq (strcat var1-label "............") 0 12))
    (when var2-label
          (setf var2-label (subseq (strcat var2-label "............") 0 12)))
    (display-string (format nil "~%~a" var1-label) w)
    (display-string (format nil ". N = ~7d Mean =~9,3f StDev =~9,3f Var =~9,3f" 
                            n1 mean1 stdev1 (^ stdev1 2)) w)
    (when var2 
          (display-string (format nil "~%~a" var2-label) w)
          (display-string (format nil ". N = ~7d Mean =~9,3f StDev =~9,3f Var =~9,3f" 
                            n2 mean2 stdev2 (^ stdev2 2)) w))
    (when (send self :paired)
          (display-string (format nil "~%Sample Differences....... Mean =~9,3f StDev =~9,3f StErr=~9,3f" meandif stdevdif sterrdif) w))
    (when (> n-ways 0)
          (display-string 
           (format nil "~%Sample Differences... DiffMean =~9,3f StErr =~9,3f Var =~9,3f (Pooled)" 
                   (- mean1 mean2) sterrdif pooled-variance ) w)
          )
   ;eliminate "exact df not known"
   ; (display-string (format nil "~2%SIGNIFICANCE TEST ~a~%Test Result."
   ;                         (if (and (not sigma) (> n-ways 0) (/= n1 n2))
   ;                             "(exact df not known)" "")) w)
    (display-string (format nil "~2%SIGNIFICANCE TEST~%Test Result.") w)
    (display-string (format nil ". ~a = ~7,3f " zt t-val) w)
    
    (when (not sigma) 
          (display-string (format nil "For df =~6,1f " (round df)) w))
    (if (< p-lev .0001) 
        (setf p-lev-string "p < .0001") 
        (setf p-lev-string (format nil "p = ~6,4f" (fuzz p-lev 4))))
    (display-string p-lev-string w)
    (when (and (not sigma) (> n-ways 0) (/= n1 n2))
   ;eliminate "exact df not known"
   ;       (display-string (format nil " (Lower Bound)") w)
   ;       (display-string
   ;        (format nil "~%                              df =~6,1f p = ~6,4f (Best Estimate)" 
   ;                (second (send self :t-dfi)) 
   ;                (send self :compute-t-or-z-p-level (second (send self :t-dfi))) ) w)
   ;       (display-string
   ;        (format nil "~%                              df =~6,1f p = ~6,4f (Upper Bound)" 
   ;                (first (send self :t-dfi)) 
   ;                (send self :compute-t-or-z-p-level (first (send self :t-dfi))) ) w)
   ;       (setf dflist (list n1 n2))
          (if (> stdev2 stdev1)
              (setf dflist (list n2 n1))
              (setf dflist (list n1 n2)))
          (setf F-test (/ (^ (max stdev1 stdev2) 2) 
                          (^ (min stdev1 stdev2) 2)))
          (setf F-df (- dflist 1)) ;(- (list n1 n2) 1)
          (setf F-p  (- 1 (f-cdf F-test (first dflist) (second dflist))))
          (display-string (format nil "~2%TEST OF HOMOGENEITY OF VARIANCE: ~%Test Result:  F = ~7,3f   df = ~s p =~9,4f" F-test F-df F-p) w))

    (unless (or sigma (not var2-label))
            
            (setf real-resp-var (if (send self :paired) "sample differences"
                                    (first real-varnames)))
            (setf effect-string (if (send self :paired) "samples" 
                                    (second real-varnames)))            
            (display-string
             (format nil "~2%STRENGTH OF RELATIONSHIP: ~%Proportion of variation in ~a related to ~a = ~4,2f" 
                     real-resp-var effect-string (/ (^ t-val 2) (+ (^ t-val 2) df))) w)  
            )

    
   ;eliminate "exact df not known"
   ; (display-string
   ;  (format nil "~2%CONFIDENCE INTERVAL: ~a~%The ~2d% confidence interval is: ~s" 
   ;          (if (and (not sigma) (> n-ways 0) (/= n1 n2))
   ;              "(exact df not known)" "")
   ;          (round (* 100 ci-level)) 
   ;          (list (fuzz (first ci-limits) 4) (fuzz (second ci-limits) 4))) w)
    (when (not sigma)
          (display-string
           (format nil "~2%CONFIDENCE INTERVAL:~%The ~2d% confidence interval is: ~s" 
                   (round (* 100 ci-level)) 
                   (list (fuzz (first ci-limits) 4) (fuzz (second ci-limits) 4))) w))
   ;eliminate "exact df not known"
   ; (when (and (not sigma) (> n-ways 0) (/= n1 n2))
   ;       (let* ((middf (second (send self :t-dfi)))
   ;              (midci (send self :compute-t-or-z-ci-limits middf))
   ;              (lowdf (first (send self :t-dfi)))
   ;              (lowci (send self :compute-t-or-z-ci-limits lowdf))
   ;              )
   ;         (display-string (format nil " for df =~6,1f" (third (send self :t-dfi))) w)
   ;         (display-string
   ;          (format nil "~%                                ~s for df =~6,1f" 
   ;                  (list (fuzz (first midci) 4) (fuzz (second midci) 4)) middf) w)
   ;         (display-string
   ;          (format nil "~%                                ~s for df =~6,1f" 
   ;                  (list (fuzz (first lowci) 4) (fuzz (second lowci) 4)) lowdf) w)
   ;         ))
    (when (send self :nonpar)
          (display-string
           (format nil "~2%~a ~a:~%Null Hypothesis:        Both Populations Identically Distributed " nonpar-st paired-st) w)
          (display-string
           (format nil "~%Alternative Hypothesis: Both Populations Differently Distributed (two tailed)") w)
          (when (> n-ways 0) (display-string
                (format nil "~%Mann-Whitney Sum of Ranks for each group = ~8,3f" 
                        (send self :mann-whitney-sum)) w))
          (display-string
           (format nil "~%~a =~5d~%Approximate  Z = ~8,3f   p = ~6,4f~%"
                   nonpar-au (send self :nonpar) (send self :nonpar-z)
                   (send self :nonpar-p)) w))
    (send w :fit-window-to-text))
  t)

(defmeth univar-model-object-proto :visualize ()
  (if (send self :var2) 
      (send (send self :anova-model) :visualize t)
     ;(send self :visualize-model-2)
      (send self :visualize-model-1))
  (setcm self)
  t)

;one sample - three side-by-side plots
 
(defmeth univar-model-object-proto :visualize-model-1 ()
  (if (not (eq current-object self)) (setcm self))
  (let* ((var (send self :var))
         (var-label (send self :var1-label))
         (obs-label (send self :labels))
         (sp (spread-plot (matrix '(2 3) 
               (list
                (histofreq var :variable-labels var-label :show nil)
                (boxplot var :diamonds t :variable-labels (list var-label)
                         :point-labels obs-label :show nil)
                (name-list obs-label :title "Observations" :show nil)
                (quantile-plot var :reg-line t :variable-label var-label
                          :point-labels obs-label :show nil)
                nil
                nil
                ))
                          :rel-widths (list 1 .8 .4)
                          :span-down (matrix '(2 3) '(1 2 2 1 0 0))
                          :model self))
         (bp (aref (send sp :plot-matrix) 0 1))
         (qp (aref (send sp :plot-matrix) 1 0))
         (hf (aref (send sp :plot-matrix) 0 0))
         (ol (aref (send sp :plot-matrix) 0 2))
         )
    (defmeth sp :spreadplot-help ()
      (plot-help-window (strcat "SpreadPlot Help"))
      (paste-plot-help (format nil "This is the SpreadPlot for single sample Univariate Analysis. Only the BoxPlot and Observation Lable windows are linked.~2%"))
      (show-plot-help)
      (send spreadplot-proto :spreadplot-help :skip t :flush nil))
    (send qp :plot-buttons :mouse-mode nil :new-y nil)
    (send bp :showing-labels t)
    (send bp :linked t)
    (send ol :linked t)
    (send bp :mouse-mode 'brushing)
    (send qp :redraw)
    (send bp :redraw)
    (send hf :redraw)
    (send ol :redraw)
    (send sp :show-spreadplot)
    sp))

;two samples - six plots
 
(defmeth univar-model-object-proto :visualize-model-2 () 
  (if (not (eq current-object self)) (setcm self))
  (let* ((var1 (send self :var1))
         (var2 (send self :var2))
         (n1 (length var1))
         (n2 (length var2))
         (var1-label (send self :var1-label))
         (var2-label (send self :var2-label))
         (paired (when (send self :paired)))
         (point-labels (combine (send self :labels) (send self :labels)))
         (sp (spread-plot 
              (matrix  '(2 3) (list
                (boxplot (list var1 var2) 
                         :variable-labels (list var1-label var2-label)
                         :point-labels point-labels :mean-line t :diamonds t
                         :connect-points paired :boxes nil
                         :show nil)
                (quantile-plot var1 :reg-line t :variable-label var1-label
                         :show nil)
               ; (histogram (combine var1 var2)
               ;          ;:variable-labels (list var1-label var2-label)
               ;          :show nil)
                (histofreq (list var1 var2)
                           :variable-labels (list var1-label var2-label)
                           :new-x t :show nil)
                (quantile-quantile-plot var1 var2 :reg-line t 
                         :variable-labels (list var1-label var2-label)
                         :show nil)
                (quantile-plot var2 :reg-line t :variable-label var2-label
                         :show nil)
                (send (send self :data-object) :obs-namelist :show nil)))
              :model self))
         (bp  (aref (send sp :plot-matrix) 0 0))
         (qp1 (aref (send sp :plot-matrix) 0 1))
         (hg  (aref (send sp :plot-matrix) 0 2))
         (qqp (aref (send sp :plot-matrix) 1 0))
         (qp2 (aref (send sp :plot-matrix) 1 1))
         (ol  (aref (send sp :plot-matrix) 1 2))
         )

    (defmeth sp :spreadplot-help ()
      (plot-help-window (strcat "SpreadPlot Help"))
      (paste-plot-help (format nil "This is the SpreadPlot for two sample Univariate Analysis. The windows in this SpreadPlot are not linked.~2%"))
      (show-plot-help)
      (send spreadplot-proto :spreadplot-help :skip t :flush nil))

    (send qp1 :plot-buttons :mouse-mode nil :new-y nil)
    (send qp2 :plot-buttons :mouse-mode nil :new-y nil)
    (send qqp :plot-buttons :mouse-mode nil :new-x nil :new-y nil)
 ;  (send hg  :plot-buttons :mouse-mode nil :new-y nil :bins t :density t)
 ;  (send hg :linked t)
    (send ol :linked t)
    (send bp :linked t)
    (send bp :showing-labels t)
    (send bp :mouse-mode 'brushing)
 ;  (send hg :mouse-mode 'brushing)
 ;  (send hg :use-color t)
 ;  (send hg :point-color (iseq n1) 'red)
 ;  (send hg :point-color (iseq n1 (+ n1 n2)) 'blue)
    (send bp :point-color (iseq n1) 'red)
    (send ol :use-color t)
    (send ol :point-color (iseq n1) 'red)
    (send ol :point-color (iseq n1 (+ n1 n2)) 'blue)
    (when (send self :paired)
          (send bp :connect-points t))
    (send qp1 :redraw)
    (send qp2 :redraw)
    (send qqp :redraw)
    (send bp :redraw)
    (send hg :redraw)
    (send ol :redraw)
    (send sp :show-spreadplot)
    sp))

;Analysis Methods

(defmeth univar-model-object-proto :analysis ()
"Args: none
Does analysis for the univariate analysis model."
  (if (= (send self :n-ways) 0)
      (send self :dependent-analysis)
      (send self :independent-analysis)))

(defmeth univar-model-object-proto :dependent-analysis ()
"Args: none
Does analysis for the univariate analysis model. Computes, for single sample or paired sample situations, the t-test or z-test value, p-level and confidence limits. When appropriate, also computes wilcoxon test."
  (send self :t-df (- (length (send self :var)) 1))
  (send self :t-value (send self :compute-t-or-z-value))
  (send self :t-p-level (send self :compute-t-or-z-p-level))
  (send self :t-ci-limits (send self :compute-t-or-z-ci-limits)) 
  (when (and (>= (send self :t-df) 7) (send self :paired)) 
        (send self :compute-wilcoxon))
  )

(defmeth univar-model-object-proto :independent-analysis ()
"Args: none
Does analysis for the univariate analysis model. Computes, for two independent sample situation, the t-test or z-test value, p-level and confidence limits. When appropriate, also computes Mann-Whitney test." 
  (send self :t-df  (- (sum (mapcar #'length (send self :var))) 2))
  (send self :t-dfi 
        (list (- (min (mapcar #'length (send self :var))) 1)
              (send self :compute-independent-df)
              (send self :t-df)))
;first df calculated as in Gravetter&Wallnau p 294---way to optimistic
;second df calculated as best estimate (Moore & McCabe p 538 (Ed 2)
;third df calculated as most conservative (Moore&McCabe p 532)
  (send self :t-value (send self :compute-independent-t))
  (send self :t-p-level (send self :compute-t-or-z-p-level))
  (send self :t-ci-limits (send self :compute-t-or-z-ci-limits)) 
  (when (>= (min (mapcar #'length (send self :var))) 5) ;7
        (send self :compute-mann-whitney))
  )

(defmeth univar-model-object-proto :compute-t-or-z-value ()
"Args: none
Computes t-test for one sample or two paired samples"
  (let ((var (send self :var))
        (mu  (send self :mu))
        (s   (send self :sigma)))
    (when (not s) (setf s (standard-deviation var)))
    (/ (- (mean var) mu) (/ s (sqrt (length var))))))

(defmeth univar-model-object-proto :compute-independent-t ()
"Args: none
Computes t-test for independent samples"
  (let* ((var (send self :var))
         (mu  (send self :mu))
         (mean-dif (apply #'- (mapcar #'mean var)))
         (nobs (mapcar #'length var))
         (pooled-var (/ (sum (mapcar #'(lambda (v) (ssq (- v (mean v)))) var))
                        (send self :t-df)))
         (std-err (sqrt (+ (/ pooled-var (first nobs)) (/ pooled-var (second nobs)))))
         )
    (send self :pooled-variance pooled-var)
    (send self :std-err-diff std-err)
    (/ mean-dif std-err)))

(defmeth univar-model-object-proto :compute-independent-df ()
"Args: none
Computes df estimate for independent samples"
  (let* ((var (send self :var))
         (variances (^ (mapcar #'standard-deviation var) 2))
         (ns (mapcar #'length var))
         (numerator (^ (sum (/ variances ns)) 2))
         (denominator (sum (mapcar (lambda (n variance) (/ (^ (/ variance n) 2) (1- n)))
                                   ns variances))))
    (/ numerator denominator)))
    


(defmeth univar-model-object-proto :compute-t-or-z-p-level (&optional t-df)
"Args: &optional t-df
Computes significance of t or z statistic for t-df degrees of freedom. Uses stored degrees of freedom if t-df nil."
  (let ((t-value (send self :t-value))
        (df (if t-df t-df (send self :t-df)))
        (direction (send self :direction)))
    (when (send self :sigma) (setf df 1000))
    (cond 
      ((not direction) (* 2 (- 1 (t-cdf (abs t-value) df))))
      ((= direction 0) (* 2 (- 1 (t-cdf (abs t-value) df))))
      ((> direction 0) (- 1 (t-cdf t-value df)))
      ((< direction 0) (t-cdf t-value df)))))

(defmeth univar-model-object-proto :compute-t-or-z-ci-limits (&optional t-df)
"Args: none
Computes confidence limits based on t or z statistic for t-df degrees of freedom. Uses stored degrees of freedom if t-df nil."
  (let* ((var   (send self :var))
         (level (send self :ci-level))
         (mean  (mean var))
         (n     (length var))
         (sterr (/ (standard-deviation var) (sqrt n)))
         (df    (if t-df t-df (send self :t-df)))
         (t-crit nil)
         (width nil))
    (when (send self :sigma) (setf df 1000))
    (when (> (send self :n-ways) 0) 
          (setf mean (apply #'- (mapcar #'mean (send self :var))))
          (setf sterr (send self :std-err-diff)))
    (setf t-crit (t-quant (- 1 (/ (- 1 level) 2)) df))
    (setf width (* t-crit sterr))
    (list (- mean width) (+ mean width))))

(defmeth univar-model-object-proto :compute-wilcoxon () 
"Args: none
Computes the wilcoxon Signed-Rank test for matched samples."
  (let* ((sign-dif  (remove 0 (send self :var) :test #'=))
         (n (length sign-dif))
         (ranks (rank-with-ties (abs sign-dif)))
         (WT nil))
    (if (< (length (which (< sign-dif 0)))
           (length (which (> sign-dif 0))))
        (setf WT (sum (select ranks (which (< sign-dif 0)))))
        (setf WT (sum (select ranks (which (> sign-dif 0))))))
    (when (not wt) (setf wt 0))
    (send self :nonpar wt)
    (send self :nonpar-z
          (/ (- wt (/ (* n (1+ n)) 4)) 
             (sqrt (/ (* n (1+ n) (1+ (* 2 n))) 24))))
    (send self :nonpar-p
          (* 2 (- 1 (normal-cdf (abs (send self :nonpar-z))))))
    t))

(defmeth univar-model-object-proto :compute-mann-whitney ()
  (let* ((vars (send self :var))
         (ranks (rank-with-ties (combine vars))) ;fwy 4.30
        ;(ranks (1+ (rank-with-ties (combine vars)))) ;ok 4.27, ng 4.28-29
         (nandm (mapcar #'length vars))
         (m (first nandm))
         (n (second nandm))
         (mn1 (+ m n 1))
         (wx (sum (select ranks (iseq m))))
         (wx2 (sum (select ranks (iseq m (1- (+ m n))))))
         (u (- (+ (* m n) (/ (* m (1+ m)) 2)) wx))
         (z (abs (/ (- wx (/ (* m mn1) 2)) (sqrt (/ (* m n mn1) 12)))))
         (p (* 2 (- 1 (normal-cdf z))))
         ) 
    (when (> u (/ (* m n) 2)) (setf u (- (* m n) u))) 
    (send self :mann-whitney-sum (list wx wx2))
    (send self :nonpar u)
    (send self :nonpar-z z)
    (send self :nonpar-p p)
    t))

;(load (strcat *vista-dir-name* "anovamob"))